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In this paper we use AP = —1.772341 ± 13.153788 s between the phenomenologically 
determined orbital period Pt, of the PSR J0737-3039A/B double pulsar system and 
C ~ ) ■ the purely Keplerian period p(°) = 2tt^/ a? /G(m\ + raj) calculated with the system's 

' parameters, determined independently of the third Kepler law itself, in order to put 

constraints on some models of modified gravity (/(/?), Yukawa-like fifth force, MOND). 
The major source of error affecting AP is not the one in the phenomenologically mea- 
5"H ■ sured period (5P^ = 4 X 10 -6 s), but the systematic uncertainty <5P(°) in the computed 

Keplerian one due to the relative semimajor axis a mainly caused, in turn, by the er- 
rors in the ratio TZ of the pulsars' masses and in sini. We get \n\ < 0.8 X 10 — 26 m — 2 
for the parameter that in the f(R) framework is a measure of the non linearity of the 
theory, \a\ < 5.5 X 10~ 4 for the fifth-force strength parameter (for A ?s a = 0.006 AU). 
J_j ' The effects predicted by the strong-acceleration regime of MOND are far too small to 

be constrained with some effectiveness today and in the future as well. In view of the 
continuous timing of such an important system, it might happen that in the near future 
it will be possible to obtain somewhat tighter constraints. 

PACS: 04.50.Kd: 04.80.Cc; 97.80.-d; 97.60. Gb 
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1. Introduction 

One of the six Keplerian orbital elements in terms of which it is possible to pa- 
rameterize the orbital motion of a pulsar in a binary system is the mean anomaly 
M. defined as M. = n(t — Tq), where n is the mean motion and To is the time of 
periastron passage. The mean motion n = 2ir/P\> is inversely proportional to the 
time elapsed between two consecutive crossings of the periastron, i.e. the anoma- 
listic period P^. In Newtonian mechanics, for two point-like bodies, n reduces to 
the usual Keplerian expression = ^GAI/a 3 , where a is the semi-major axis of 
the relative orbit of the pulsar with respect to the companion and M = m p + m c 
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is the sum of their masses. In pulsar timing the period Pb is very accurately de- 
termined in a phenomenological, model-independent way, so that it accounts for all 
the dynamical features of the system, not only those coming from the Newtonian 
point-like terms, within the measurement precision. 

In this paper we wish to look for phenomenologically determined deviations 
from the third Kepler law in the PSR J0737-3039A/B system^, consisting of two 
radio-pulsars moving along a moderately eccentric, 2.4 hr orbit (see Table [T] for 
its relevant orbital parameters), in order to dynamically constrain some models 
of modified gravity {f(R), Yukawa-like fifth force, MOND). We believe that the 

Table 1. Relevant orbital parameters of the PSR J0737-3039A/B system 2 . The projected 
scmimajor axis is defined as x = (a^ c /c) sini, where a^c is the barycentric semimajor axis, 
the phenomenologically determined post-Keplerian parameter s is equal to sin i in general 
relativity, and i is the angle between the plane of the sky, perpendicular to the line-of-sight, 
and the orbital plane. We have conservatively quoted the largest error in s reported in 
Ref. 2. The parameter e is the eccentricity. The orbital period is known with a precision 
of 4 X 10~ 6 s. The quoted values of m A and mj were obtained in Ref. 2 by using the 
general relativistic expression of the post-Keplerian A's periastron advance for the sum 
of the masses M = 2.58708(16), and the ratio of the masses TZ = m A /m-g, = 1.0714(11) 
phenomenologically determined from both the projected semimajor axes ia and xb . 
The rotational period V A = 2n/n A of PSR J0737-3039A/B A amounts to 22 ms, while 
V B = 2w/n B = 2.75 s. 

Pb (d) x A (s) s e m A (M Q ) ra B (M Q ) 

0.10225156248(5) 1.415032(1) 0.99974(39) 0.0877775(9) 1.3381(7) 1.2489(7) 



analysis presented here has the merit of quantitatively determining how useful a 
laboratory such as PSR J0737-3039A/B and a method like the orbital period can 
be in testing alternative theories of gravity by precisely individuating the major 
sources of errors. Indeed, as many different laboratories and methods as possible 
should be used to extensively tests a fundamental interaction as gravity. Such an 
analysis may also yield insights concerning future obtainable improvements in view 
of the continuous timing of such a binary system. 

2. Deviations from the third Kepler law 

Looking for genuine deviations from the third Kepler law— meant as discrepancy 
AP between the phenomenologically determined orbital period Pb and the purely 
Keplerian one P<°) calculated with the values of Table [TJ-from the PSR J0737- 
3039 A/B timing data is possible because all the parameters entering AP have been 
just measured independently of the third Kepler law itself. Indeed, the relative 
semimajor axis entering P' ' 

a = (l+H) (^) (1) 

is built in terms of the ratio 1Z of the masses, the projected semimajor axis xa and 
sini; the projected semimajor axes xa and xb were phenomenologically determined 
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from the timing data, the phenomenologically estimated post-Keplerian parameter 
s can be identified with sin i in general relativity, the ratio 1Z has been phenomeno- 
logically determined from the ratio of the projected semimajor axes cor ning from 
the quite general relation, valid to at least the first post-Newtonian ordep214l 



(A) (B) 



(2) 



and the sum of the masses has been, in turn, derived from the A's periastron 
advance which is, at present, the best determined post-Keplerian parameter. The 
third Kepler law only enters the expression of the mass functions 

77ig^ A sin 3 i 

•^A/B = 7 j 72 : ( 3 ) 

which, instead, have not been used in obtaining the parameters of PSR J0737- 
3039A/B . 

Table Q] yields 

AP = P h - P (0) = -1.772341 ± 13.153788 s. 
In Eq. ([4]) we have evaluated 



S(AP) < 6P h + SP( 0) , 



with 



6P(°)<SP(% + 6PM\ U + 6PM\ G . 
The terms in Eq. ([6]) are 



(4) 
(5) 
(6) 



( 5P^\ a < 3TTy/^Sa = 12.218572 s, 



SP {0) \m <= n^Jc^SM = 0.273248 s, 
SP {0) |g <= ttJ-^SG = 0.661964 s, 



(7) 



so that 



(8) 



Sp(o) < 13.153784 s. 

In Eq. we used the values of Table Q] for SM ancP 6G = 0.0010 x lO" 11 kg" 1 
m 3 s^ 2 ; the uncertainty in a, which is the most important source of error, has been 
evaluated as 



ith 



5a < Sa\ n + Sa\ s + Sa\ XA , 

' Hn < (^) $K = 466, 758 m, 

Sa\ s <(1+TZ) (ff ) 5s = 342, 879 
1<H»a < (1+11) (f) fa A = 621 m. 



(9) 



(10) 
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Thus, 

8a < 810, 259 m. (11) 

It is interesting to note that 1Z and s have a major impact on the overall uncertainty 
in a; our estimate has to be considered as conservative because we adopted for Ss 
the largest value quoted in Ref . [5] 

The following considerations about the treatment of errors are in order. 

• Realistic upper bounds of the errors have been obtained by purposely lin- 
early summing the various sources of uncertainty because of the existing 
correlations among the estimated parameters of the PSR J0737-3039A/B 
systenP 

• The uncertainties of Table Q] and used throughout the paper are twice the 
parameter uncertainties given by the software used in Ref. [2 the authors 
of such a work believe that the real measurement uncertainties are actually 
somewhat smaller than those quoted. If so, also our result for AP would, 
in fact, be more accurate. By the way, in view of the continuous timing of 
the PSR J0737-3039A/B system it is likely that in the near future the pre- 
cision reached in determining its parameters will allow to better constrain 
AP, although it is difficult to quantify the extent of such an improvement 
because of the impact of xb via 1Z. 



3. Applications 

In this Section we will look for deviations from the third Kepler law induced by small 
radial accelerations A r , in general non-constant, induced by some modified models 
of gravity (f(R), Yukawa-like fifth force, MOND). The effects of such accelerations, 
which are to be retained quite smaller than the Newtonian monopole, should, in 
principle, be treated perturbatively with, e.g., the Gauss perturbing equations for 
the variations of the Keplerian orbital elements of a test body. Indeed, while the 
average Newtonian acceleration for the PSR J0737-3039A/B system amounts to 

CM 

(A N ) = — 7 == = 446 ms- 2 , (12) 

the f(R) acceleration (see Section [3~Tj) is about three orders of magnitude smaller, 
the Yukawa-type term (see Section I3.2[) is five orders of magnitude smaller, while 
the MONDian correction (see Section [373]) is even much smaller. 

However, in view of the likely broad spectrum of readers, many of whom may 
not be fully acquainted with the rather specialistic tools of celestial mechanics, we 
will simply work with 

<j> 2 r = — + A r , (13) 

where cj> is an azimuth angle, by neglecting the orbital eccentricity since the small 
corrections due to it would not affect the constraints set by the main terms. 
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3.1. Extended theories of gravity: "f(R)" 

It has been suggested that cosmic speed-up can be explained by generalizing general 
relativity introducing in the gravitational action terms non linear in the scalar 
curvature R (see Refs. [6] [7] and references therein). These theories are the so called 
f(R) theories of gravity (or "extended theories of gravity"), where the gravitational 
Lagrangian depends on an arbitrary analytic function / of the scalar curvature R. 
If, on the one hand, extended theories of gravity seem to be appealing in cosmology, 
on the other hand, there is the need of checking f(R) predictions with Solar System 
tests which general relativity passed with flying colors^. Actually, the debate on the 
compatibility of these theories with Solar System tests is still open (see Ref. [9] and 
Refs. [lOl [IT] and references therein). Here we do not enter the debate, but simply 
suggest how, on the basis of the approach already outlined in Refs. [TU1 [HI the data 
from the double pulsar can be used, at least in principle, to constrain the allowable 
analytical forms of the f(R). We work in the Palatini formalism (for an approach to 
extended theories of gravity in the metric formalism see, for instance, Ref. [12] and 
references therein). The field equations of f(R) theories, in the Palatini formalism, 
which have the explicit form 

f(R)R {flv) (T) - \f{R)g^ = nT™ at (14) 

V^[V5/'(W y )=0 (15) 

admit the spherically symmetric vacuum exact solution: 

, 2 / 2GM kr 2 \ 2)2 dr 2 
* =-{ l - — + —) c d * + (! _ 2gM + 

+ r 2 dd 2 + r 2 sin 2 M(p 2 , (16) 

which is referred to as Schwarzschild-de Sitter metric (see Ref. [T3] for details on the 
derivation of this solution) . The parameters appearing in Eq. (|16|) are the mass M 
of the spherically symmetric source of the gravitational field and k, which is related 
to the solutions R — c.i of the structural equation 

f(R)R-2f(R)=0. (17) 

controlling the solutions of Eq. (|14p . Namely, it is k — Cj/4 = R/A. Indeed, we may 
say that k is a measure of the non-linearity of the theory (if f(R) = R, Eq. (fT7|) 
has the solution R = — > k — 0) In general relativity, the Schwarzschild-de Sitter 
solution corresponds to a spherically symmetric solution of Einstein equations, with 
a " cosmological" term Ag M „, A being the cosmological constant. In practice, it is 
A = — k in our notation. As we can see from Eq. (fTB"]) . the modifications to the 
solutions of the field equations due to f(R) theories are given by a term proportional 
to the Ricci scalar. In what follows, we write the gravitational potential in the form 



U = U Q + U K , 



(18) 
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with 

U = -^-, 

(19) 

U K = -KT 2 

where n — fc/3, and we consider AU = U K as a perturbation. In doing so, we 
neglect the effect of spatial curvature. From the perturbing potential Uk in Eq. ([T9"]) 
we obtain an entirely radial perturbing acceleration 

A = -2/cr. (20) 

We may write 

P b = P<°) + P< fc ) (21) 

where the contribution AP = p( k \ due to the non linearity of the Lagrangian, by 
applying Eq. (fT3|) . turns out to be 

AP = -^, (22) 

from which k can be evaluated 

AP "'°' 3 (23) 



2tt 

Consequently, we may give the following estimates 

fc <^ (AP)+ 3|A^'°'V (24) 

~ 2tt v ; 2tt v ; 

and 

k= (1±7) x 10" 10 s^ 2 (25) 
Equivalently, restoring physical units 

K = (0.1 ±0.8) x 10" 26 m~ 2 , (26) 
or, rephrased in terms of the cosmological constant, 

|A| =3k= (0.3 ±2.5) x 10~ 26 m" 2 (27) 



The estimates of Eq. (1271) can be used to set constraints on the analytical form of 
the f(R), according to an approach described in Refs. QUI EJ For instance, if we 
consider the Lagrangian 

f(R) = R-£. (28) 

which mimics cosmic acceleration (even if it has well known instabilities see, e.g., 
Ref. I14p . thanks to the estimates on k we may set a limit on the parameter /i. We 
get 

\fi\ < 10~ 20 eV. (29) 
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We remark that the estimates of Eq. ([2"7) are of the same order of magnitudes of 
those obtained in Refs. [10l [Til where data coming from Solar System observations 
were used. So, as we did there, we can say that also data coming from pulsars are 
hardly fit to constrain f(R), since the value Eq. (|29| is remarkably greater than 
estimate^ /x ~ 10 -33 eV, needed for f(R) gravity to explain the acceleration of the 
Universe without requiring dark matter. 



3.2. Yukawa-like fifth force 



Modifications of the gravitational potential leading to a Yukawa-like term arise in 
different contexts, such as brane- world models, scalat-vector-tensor theories of grav- 
ity, string theory, study of cosmological defects (see, e.g., Refs. \15\ [TBI and references 
therein) . In general, a Yukawa- like force can be obtained from the potentiaff^ 



U = U + U Y , 



with 



U = - 



GM 



(30) 



(31) 



^ = -^[aex P H)] ) 

where |a| is the strength of the Yukawa- like force and A its rangJS. The parameters 
a, A have been bounded thanks to many laboratory experiments and observations, at 
different scales (see Ref. [19] and references therein) . Here we show how the strength 
\a\ of the Yukawa-like force can be constrained at a scale A « a « 10 9 m, which is the 
length scale of the double pulsar PSR J0737-3039A/B . The perturbing potential 
AU = Uy in Eq. ([3Tj) leads to the following expression for the radial component of 
the corresponding perturbing acceleration: 

GMa 



As 



( 1+ D exp H) 



From Eq. 



it is possible to obtain, in the circular orbit approximation, 



P = p(°) 



l-SL 

2 



-,exp 



Hence, on writing 



(32) 



(33) 



(34) 



from Eq. (113|) we obtain for the contribution AP = p( Y ) the following expression 



AP 



na / r 
# V + A 



exp -- 



By assuming A ; 



r w a, we get 

APn (0) cxp(-l) 
2tt 



= 0.7 X 10" 



S a < n^EUl s(AP) 



(35) 



(36) 



|AP|exp(-l) J (0) 

Z7T 



5.5 x 10" 
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Thes e estimat es for a, should be compared with those obtained at different 
scale ;p0 | 21 | 22 | j n particular, in Ref. [21] it is showed that |a| < 10~ 9 at a scale 
corresponding to that of the double pulsar. In other words, our results are not com- 
pelling, since they are orders of magnitude greater than the best estimates already 
available. This means that the length scale of the hypothesized Yukawa-like fifth 
force is much greater or much smaller than the length-scale of the double pulsar 
PSR J0737-3039A/B (a = 0.006 AU). The results obtained in Ref.E2for the Solar 
System point towards the second hypothesis, i.e. A < 0.006 AU. 

3.3. MOND 

MOND predicts that the gravitational acceleration A g felt by a particle in the field 
of a distribution of mass is 

A 9 = (37) 

where An is the Newtonian acceleration, Aq is a n acceleration scale which different, 
independent ensembles of observations set tcE^ An — 1.2 x 10~ 10 m s -2 , and fi(x) 
is an interpolating function which approximates 1 for x 3> 1, i.e. for accelerations 
larger than An; for x <C 1 = x, so that in such a strongly MONDian regime 
A g ps \/ An An. For a quite general class of interpolating functions, jJi{x) can be cast 
into the formal 

, (38) 
which yields a modified gravitational acceleration^! 



A g « A N 1 + fc 
Note that the most commonly used expressions for f/,(x), i.e. 



(39) 



n(x) = -=4=, (40) 
V 1 + x 



= —-, (41) 
1 + x 

can be obtained from Eq. (|3"9")) for fc = 1/2, q — 2 and k = 1, q = 1, respectively. 
A useful review dea ling, among other things, with many attempts to theoretically 
justify MOND i gain, we may write 

P b = p(°)+P( M \ (42) 

and, from Eq. (fT3"|) with A r given by Eq. (|3"9"|) , on setting e — ► 0, it is straightforward 
to obtain for the perturbation AP = p( M \ 



P^ko (A, 



AP=- ; -(-^) . (43) 
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By defining 

it turns out that both Eq. (gDJ) and Eq. (gTJ) are compatible with the PSR J0737- 
3039 A/B data: indeed, 

■Hi/2 2 = (-1.610 ±0.008) x 10~ 22 s, 

(45) 

= (-1.431 ± 0.004) x 10~ 19 s. 

By the way, from Eq. ([15]) it is clear that the PSR J0737-3039A/B system is defi- 
nitely unsuitable for testing MOND because discrepancies as small as those of Eq. 
(|45p will never be measurable, falling well below even the precision in determining 
Pb- The same holds also for other binary systems, given the obtainable accuracy in 
determining their orbital periods. Solar System dynamicJ^ has proven to be more 
effective in constraining the interpolating function for x » 1. 



4. Discussion and conclusions 

In this paper we looked for deviations from the third Kepler law in the PSR J0737- 
3039A/B double pulsar system finding a discrepancy AP = -1.772341 ± 13.153788 
s between the phenomenologically determined orbital period Pb and the purely 
Keplerian one p(°\ While Pb is determined to a 10 -6 s level, P^ can be considered 
known only at w 10 1 s level after a conservative propagation of the uncertainties 
in the system's parameters entering its expression. The consistency of our analysis 
is assured by the fact that the values of the orbital parameters involved in the 
calculation of P" have been determined independently of the third Kepler law 
itself. The major source of error in p(°) is due to the semimajor axis which, in turn, 
is mainly affected by the ratio 1Z of the pulsars' masses, which, in turn, depends 
on the poorly determined projected semimajor axis of B, and by sini. Continuous 
timing of the PSR J0737-3039A/B system might reduce such errors in the near 
future, perhaps yielding more precise results. 

We used AP to constrain various theories of modified gravity like f{R), Yukawa- 
type interactions and MOND. The bound on the parameter k, that in f(R) 
framework is a measure of the deviation of the theory from general relativity, is 
\k\ < 0.8 x 10~ 26 m -2 ; differently speaking, the bound on k can be seen as a bound 
on the cosmological constant |A| = 3k. By assuming a range A of the order of 
the PSR J0737-3039A/B system's size, we get |a| < 5 x 10~ 4 for the strength of a 
putative Yukawa-like fifth force. In regard to MOND, the PSR J0737-3039A/B sys- 
tem is neatly unsuitable to test it because the order of magnitude of the predicted 
MONDian discrepancies is AP M qnd ~ 10~ 19 - 10~ 22 s. 
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